% Computes bi-variate transfer functions between two channels, given
% two input & one output channel
% Latest date 21.11.2002
% bug fixed - correct signs for the cross spectra
% compared with readedi1 erros, z coh found correct

%zij = za zb - zc zd;
%tf = zxy zyx zxx zyy

function[tf] = jktf(SPMall,ProcDef,iflag),

SPM(:,:,:) = SPMall(1,:,:,:);

nfreq=ProcDef.nfrq;
EQN=ProcDef.equations;
neq=ProcDef.neq;

tf = zeros([ProcDef.len neq]);

for l = 1:neq,
i = EQN(l,1);
j = EQN(l,2);
k = EQN(l,3);
p = sqrt(-1);

za = SPM(:,j,j);
if i>k,
   zb = SPM(:,i,k)+p*SPM(:,k,i);
elseif k>i,
   zb = SPM(:,k,i)-p*SPM(:,i,k);
end;

if k>j,
   zc = SPM(:,k,j)-p*SPM(:,j,k);
elseif k<j,
   zc = SPM(:,j,k)+p*SPM(:,k,j);
end;

if i>j,
   zd = SPM(:,i,j)+p*SPM(:,j,i);
elseif i<j,
   zd = SPM(:,j,i)-p*SPM(:,i,j);
end;

den = SPM(:,j,j).*SPM(:,k,k)- (SPM(:,k,j) - p*SPM(:,j,k)).*(SPM(:,k,j) + p*SPM(:,j,k));
top = [za.*zb-zc.*zd SPM(:,k,k).*zd-conj(zc).*zb];

tf(:,l) = (top(:,1)./den).*iflag;
tf(:,l+2) = (top(:,2)./den).*iflag;

end;